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Abstract. A model for polar filaments interacting via molecular motor complexes is investigated which 
exhibits bifurcations to spatial patterns. It is shown that the homogeneous distribution of filaments, such 
as actin or microtubules, may become either unstable with respect to an orientational instability of a finite 
wave number or with respect to modulations of the filament density, where long wavelength modes are 
amplified as well. Above threshold nonlinear interactions select either stripe patterns or periodic asters. 
The existence and stability ranges of each pattern close to threshold are predicted in terms of a weakly 
nonlinear perturbation analysis, which is confirmed by numerical simulations of the basic model equations. 
The two relevant parameters determining the bifurcation scenario of the model can be related to the 
concentrations of the active molecular motors and of the filaments respectively, which both could be easily 
regulated by the cell. 



PACS. 87.16.-b Subcellular structure - 47.54.+r Pattern formation 



.75.-k Complex systems 



1 Introduction 

In eukaryotic cells the polar filaments actin and micro- 
tubules interacting with motor proteins play a crucial role 
in intracellular organisation and transport as well as for 
the static and dynamical structure of the cytoskeleton ^ 
12]. Most prominently, microtubules and kinesin are in- 
volved in highly connected dynamical structures, such as 
the mitotic spindle in cell division [30], while the motility 
of the cell as a whole is governed by acto-myosin com- 
plexes . The dynamical behavior of such assemblies 
depends on the available biological fuel ATP (Adenosine- 
triphosphate), which is consumed during the motion of 
motor proteins along the filaments. Vesicles for instance 
are transported across a cell by motors moving along the 
tracks defined by microtubules, or - which is the scope of 
this work - oligomeric motor proteins that attach to two 
or more filaments induce relative motion between neigh- 
boring filaments and cause dynamical networks. The latter 
process is vitally important in cells, since the cytoskeleton 
constituted of the filaments has to be self-organized and 
even actively reorganized not only during cell-locomotion 
but also in order to react to outer stimuli (chemotaxis). 
During mitosis, microtubules attach to the chromosomes, 
which are then divided and the two halves finally are 
transported by the motor-induced filament sliding into 
the two evolving daughter cells. 

Since the situation is very complex in a living cell, 
well designed in vitro experiments are the agent of choice 
for controlled explorations of prominent aspects of cellu- 
lar systems. Recent experimental progress yielded indeed 



important insights into organization and dynamical prop- 
erties of the cell, which in turn call for modeling activities 
to foster their deeper understanding. These experiments 
comprise investigations of self-organization in filament- 
motor mixtures of microtubules in the presence of a single 
type of motor protein 0IHE] and more recently also in 
actin-myosin networks |1(JI11| as well as assays where two 
types of motors interact with microtubules J2]- Even in 
such model systems, simple compared to a living cell, there 
has been found a great variety of different two-dimensional 
patterns, such as stripe patterns, asters, vortices and ir- 
regular arrangements. 

Works on modeling active filament motor systems com- 
prise molecular dynamics simulations |51ll2U13| and mean 
field approximations for spin-like models displaying stripe- 
like patterns ^1] . More recently a phenomenological model 
of a motor density interacting with a phenomenological 
vector field has been considered which is able to reproduce 
asters and vortex-like solutions as wen as models for 
bundle contraction in one spatial dimension |16II17II18| . 

Here we follow a more microscopic approach starting 
from the Smoluchowski equation for the spatial and an- 
gular distribution of rigid rods which approximate 
the rather stiff microtubules and with limitations also the 
actin filaments. It is well known since Onsager's theory 
that an ideal rigid rod system exhibits a nematic liquid 
crystalline order beyond some critical rod density |2l)U21) 
which is still present in generalizations for semiflexible 
polymers, albeit at a slightly higher filament density \'22\. 
and has been observed in vitro for both kinds of biopoly- 
mers |23U24| . To use Onsager's equilibrium argument of 
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excluded volume induced transitions in a nonequilibrium 
system like the cell, one can formulate phenomenologi- 
cal models as done recently to describe how nucleation 
and decay of the filaments in a cell at a density close 
to the isotropic-nematic transition may give rise to spa- 
tially periodic patterns . Alternatively one can use the 
aforementioned Smoluchowski equation for rigid rods 
which is an approach widely used in polymer and colloidal 
science and recently has been supplemented by active cur- 
rents to describe filament -motor-systems |26) in a way 
inspired by a macroscopic model for filament bundling in 
one spatial dimension |17| . In this approach one makes 
use of the fact that the small motor proteins diffuse much 
faster than the filaments. Hence the density of the mo- 
tors can be assumed to be homogeneous and - as well 
as properties like the mean velocity and the duty ratio 
of the motors - enters into the model only via the coef- 
ficients. The phcnomcnological description of the active 
currents can be derived by symmetry considerations and 
includes three major contributions, whereof in this work 
we focus on the effects on pattern formation of only two 
of them. Since most of the experimental assays are quasi 
two-dimensional, we restrict our analysis also to two spa- 
tial dimensions. 

The filament-filament interactions induced by the mo- 
tors as well as by excluded volume effects are nonlocal. 
Thus they have to be approximated by a gradient expan- 
sion, which must not be truncated at the leading order, 
because on their basis the instabilities of the homogeneous 
filament distributions cannot be predicted as pointed out 
recently [23, but has to be continued up to fourth or- 
der derivatives. Since in vitro aster-like patterns evolve at 
much lower filament density than the isotropic-nematic 
transition, a moment expansion of the probability density 
function can be truncated to derive a closed set of equa- 
tions for the physical observables, which are the density 
and the orientation field of the filaments. After a detailed 
linear analysis the major goal of this work is a character- 
ization of the nonlinear behavior of patterns beyond their 
threshold. For this purpose we employ on the one hand 
numerical simulations of the coupled equations for the fil- 
ament density and the orientational field. On the other 
hand we use the powerful method of amplitude expan- 
sion, where equations of motion for the amplitudes of the 
spatially periodic pattern are derived close to the pattern 
forming instability [5%ll29ll3()| . 

The generic patterns in two-dimensional systems close 
to threshold are stripes, squares or hexagonal patterns and 
for the amplitudes of these patterns generic nonlinear am- 
plitude equations may be derived |28II31| . In numerical 
simulations we found either stripes or squares (the absence 
of hexagonal structures can be explained by the hierarchy 
of equations in the derivation of the amplitude equation) 
and therefore we focus on the equations for stripe and 
square patterns (as described by Eq. (|?7jl below), which 
are of the form |2"£ll3"2"] 

T d t X = eX- 9l \X\ 2 X - g 2 \Y\ 2 X , (la) 
r d t Y = eY - 9l \Y\ 2 Y - g 2 \X\ 2 Y . (lb) 



Here e denotes the deviation of the control parameter from 
its value at the threshold of pattern formation, while X 
and Y describe the amplitudes of a spatially periodic pat- 
tern in one (stripe) or two orthogonal directions (square), 
namely the x- and y-direction. These amplitude equations 
are universal for all pattern forming systems of a given 
symmetry and the specific information of each system is 
only encoded in the coefficients, namely To, 51,(72, which 
are functions of the parameters of the underlying system 
under consideration. The nonlinear analysis of Eqs. Q 
is worked out analytically and by this advantage a whole 
phase diagram for the various patterns close to threshold is 
presented. In contrast, such an investigation would be very 
cumbersome and time consumptive using numerical sim- 
ulations of the model equations. The solution \X\ = \Y\ 
corresponds to a square pattern for the orientational field 
of the filaments which resembles very much a structure of 
periodically arranged asters. 

The work is organized as follows. The underlying mi- 
croscopic model for the filament distribution is presented 
in Sec. El and in the same section we derive also the cou- 
pled equations for the density and orientation of the fila- 
ments, whereby the technical parts are moved to the ap- 
pendices. The thresholds of the pattern forming instabili- 
ties and its dependence of the parameters are determined 
in section |3 where we find an instability leading to den- 
sity modulations, an orientational instability as well as a 
liquid-crystalline isotropic-nematic (I-N) transition. With 
regard to the formation of asters, the orientational insta- 
bility is the most interesting one and its nonlinear behav- 
ior, namely the nonlinear stripe solutions and the periodic 
lattice of asters and inverse asters as well as their stability 
are investigated in more detail in Sees. 0] and The work 
is finished with a summary and concluding remarks. 

2 The model 

The probability of finding a rod (with fixed length L) at 
the position r with the orientation u (with |u| = 1) at 
time t is described by the distribution function </7(r, u, i) 1 
and the governing so-called Smoluchowski equation is just 
the continuity equation for the probability 

dtW + V • J t + K ■ J r = . (2) 

Homogeneously distributed molecular motors are supposed 
to interact with the rods inducing active currents on them. 
The total translational and rotational currents are given 
in Cartesian coordinates by 

Jt,i = -Da [dp + SPfyVex] + J? ti , (3a) 
J r<i = -D r [Hi* + <PK t V ex } + J«i . (3b) 

The terms with an upper index a are the motor-mediated 
active currents, the translational diffusion matrix reads 

Dij = D\\UiUj + Z)_l (Sij - UiUj) , (4) 

1 In the following we will write ^(r, u) for reasons of brevity 
and ^(r, u, t) only if we want to emphasize the time depen- 
dence. 
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D r is the rotational diffusion constant and the operator of 
rotational diffusion is given by 



K = u x 3 U 



V PT describes the excluded volume interaction 



(5) 



14x(r,u) = J du 1 J dr 1 VK(r-r',u,u'Mr',u') (6) 

between rods, where the interaction kernel is expressed in 
terms of the so-called Straley coordinates £ and r\ in two 
dimensions 

L/2 L/2 

W(r-r',u,u') = |uxu'| J dC, J drj S (r-r'+uC+u'^) . (7) 

-L/2 -L/2 

This expression takes into account that there is only an in- 
teraction between rods with coordinates (r, u) and (r', u') 
in the case of a finite overlap i.e. if the connection vector 
r — r' can be constructed by a linear combination u^ + u'rj 
of the rod orientations with —L/2 < £, -q < L/2. The ac- 
tive currents induced by the homogeneous motor density 
are 

Jt = <?(r, u)J du'J dr'v(r-r' u, u')W(r-r', u, u')W(r', u') , 
J r a = <?(r, u)Jdu'J dr'u(u, u')W(r-r', u, u')&(r', u') (8) 



with the translational and rotational velocity given by 



, . a r — r' 1 + u • u' (3 u' u 
v(r-r,u,u ) = -- 



2 L luxu' 



2 |u x u' 



w(u,u') = 7(u • u'). 



U X u' 



(9) 



(10) 



respectively. The active currents are normalized to the in- 
teraction volume and they fulfill both, the conservation of 
translational and rotational momentum in the absence of 
external forces and torques, as well as translational and 
rotational invariance (c.f. |2H] and appendix IX}) . 

According to its functional dependence 1 + u • u', the 
first term in Eq. @ contributes at maximum to the fila- 
ment transport for a parallel orientation, which identifies 
it as the leading pattern forming mechanism in counterac- 
tion with the diffusive motion. It arises from a difference 
in motor activity between the filament center and the fila- 
ment endpoints, where motors stall for some finite time 77, 
which has been recognised as a crucial condition for aster 
formation |H4| . It depends also on the relative filament 
separation r — r', hence no relative motion is obtained if 
the motor connects two filaments having the same centers 
of mass. 

The second term in Eq. @, proportional to (3, causes 
filament translations along u' u and therefore becomes 
maximal for antiparallel filament orientations. The contri- 
bution in Eq. I|1U|I causes rotational currents. In this work 
we focus mainly on active currents caused by the a-term 



and especially on its effects on pattern formation in two 
spatial dimensions. 

Eq. (J2J together with the currents as defined by Eqs. 10 
and (JSJ) cannot be solved analytically for ^(r, u, t) and 
even numerically, being a nonlinear integro-differential e- 
quation, it would be rather intricate. However, these equa- 
tions may be simplified by expressing the distribution func- 
tion in terms of its moments. For a rod distribution, the 
zeroth, first and second moment with respect to the orien- 
tation correspond to the filament density p(r,t), the po- 
lar orientation t(r, t) and the nematic order parameter 
Sij (r, t) , which are defined by the following expressions 

p(r,t) = J du &{r,u,t) , 
t(r, t) — J du u \P(r, u,t) , 
Sij(r,t) = J du Ul u 3 ;#(r,u,i) . (11) 

In contrast to a usual lyotropic liquid crystal [201 1 here the 
filaments are polar with respect to the motor action which 
breaks the ±u-symmetry and therefore the first moment 
may become finite. 

The critical filament density of the homogeneous iso- 
tropic-nematic (I-N) transition can be determined by the 
stability of the second moment of the rotational diffu- 
sion contribution, i.e. where — D r J u a upTZ [TZ'P + \I/lZV ex \ 
changes its sign. This condition yields the critical density 



Pin = 2 71 " 



(12) 



in two spatial dimensions where the I-N transition is of 
second order. 

In order to obtain a closed set of equations of motion, 
the distribution function ^(r, u, t) may be expanded with 
respect to the leading moments as given by Eqs. 1)11(1 . 
In the parameter range well below the homogeneous I— 
N transition, where aster formation is observed in vitro 
and on which we concentrate in this work, the first two 
moments are sufficient and the distribution function may 
be represented in terms of pit, t) and t(r, t) (for more 
details we refer to appendix IO) 

!P(r, u, t) = ^ (p(r, t) + 2u ■ t(r, t)) . (13) 

To get rid of the integral kernel in Eq. Q, a gradient 
expansion can be performed as described in more detail 
in appendix [5] After a rescaling of variables via 



Du 1 . 

t = ~p t ' x = J x ' p = VqP ' * = Vot 

I 



, I 2 , D ± 1 

a ' = tt" ' D r = TT r ' d = TT = o 

L)\\ JJu U\\ I 



(14) 



where for the last relation we refer to ^1, we obtain in two 
spatial dimensions the following coupled equations for the 
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density and the two components of the polar orientation 
field ti (i — x,y) of the filaments 



therefore influences only the nonlinear behavior of pattern 
selection beyond threshold as shown briefly in Sec. 0] 



dtp = ^Ap 



-rrft 



Ci 
16ft 



3 «]_._. 
2tt 24j yl y ' 



38V • (pVAp) + lift (tjdiAtj 



tiAdjtj + 2tjdjdiditi + tjdjAti 



, (15a) 



5 1 
d t U = -D r U + -AU + 7ft v • t 
8 4 



1 

4tt 

-ft 

96 3 

— ft 
96 



5ft (tidjp) + djitjdip) + difodjp) 



3tidjp + tjdip + p{ditj + djU) 
16a 



tidip + pditi 



C 2 



-ft 



pAdiU + tidiAp 



p[ lid j At i + 16diAtj + 32ft ft ft i ; 



16t j d l Ap + 32t ; ftftftp + UtidjAp 



1 

48 



2 -Id 

4 7T ? 



tjdjdip - -UAp 



(15b) 



with Ci = 23040 and C* 2 = 2Ci. 

The contributions to Eqs. i|15[l that include also fourth 
order derivatives are crucial for the determination of the 
instabilities evolving from the homogeneous filament dis- 
tribution and the nonlinear patterns in the following sec- 
tions. Their importance was pointed out recently in a com- 
ment on a previous work [23 . It should be noted that 
fourth order derivatives proportional to the translational 
diffusion have been neglected here because they are over- 
compensated by the fourth derivatives resulting from the 
active current and therefore influence the presented results 
only quantitatively by a small amount (but considerably 
complicate the equations). 

The instabilities of the homogeneous filament distri- 
bution caused by the a-contribution to the active current, 
c.f. (|9I10|) . and the nonlinear solutions are determined in 
the following sections. The active current proportional to 
(3 leads to additional effects, which are responsible e.g. for 
propagating modes in the system and for the breaking of 
the unexpected t — ► — t symmetry of Eqs. (|15(l . To remind 
the reader, we have allowed t to become nonzero, so we al- 
lowed the breaking of the ±u-symmetry of Eq. (J2J, but our 
model equations (|15|l nevertheless have a it-symmetry, 
since terms like t 2 appear in the equation for t only due to 
the /3-contribution of the current. These effects are deter- 
mined elsewhere |35| and the corresponding contributions 
to Eqs. H15|) are not shown for the sake of brevity. The ac- 
tive rotational current proportional to 7 is nonlinear and 



3 Threshold for density and orientational 
instabilities 

Having defined the model, the first issue in the field of 
pattern formation is the determination of possible insta- 
bilities. As calculated in this section by a linear stability 
analysis, in a certain parameter range the homogeneous 
basic state, which consists here in a constant filament 
density po and a vanishing polar orientation t = 0, be- 
comes unstable with respect to inhomogeneous perturba- 
tions p(r,t) and t(r,i). For this purpose, by the ansatz 
p(r, t) = po + p(r, t) we separate the constant part po 
of the filament density from the spatially inhomogeneous 
one, p(r,t), and linearize Eqs. I|15|l with respect to small 
inhomogeneous contributions p(r,t) and t(r,i). Accord- 
ingly a set of three coupled linear equations is obtained 



ftw(r, t) = £ w(r,t) 



(C^ 


V 



r (o) r (o) 

'-22 i -23 

r (0) AO) 

i -32 i -33 , 



r(v,t) , (16) 



for the three components of the vector 



ir(r,t) = t* (r,i). 



(17) 



The components of the linear operator £0 are given by 
the expressions 



,(0) 



-11 



C 

Con = —D 



3 / 2 \ ap 

4 [ n P °) ~~2A 



A 



19 apo 
11520 ' 



apa 
46080 



-(0) 
-23 



(11Z\ 2 + 64A9 2 ) 



48 J 



720 



Adrd., 



(18) 



and the two further components and £33^ may be 

obtained by permuting ft and ft, in and C 1 ^ , re- 
spectively. Naturally the mode ansatz 



w(r, t) = E exp (at + ik • r) 



(19) 



with the wave vector k = (q,p) and the eigenvector E 
solves the linear homogeneous set of equations (TIT))) and 
the solubility condition provides a third order polynomial 
for the eigenvalue a, which factorizes into a linear and a 
quadratic polynomial describing different types of insta- 
bilities. Considering moderate filament densities, for in- 
termediate values of a, an orientational instability with a 
finite wavelength is preferred, whereas the density mode 
does not couple on the level of the linear equations. For 
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large values of a a further mode becomes unstable which 
resembles a spinodal decomposition of the filament density 
driven by the motor proteins. 

The latter instability with respect to density fluctua- 
tions is governed by and the eigenvalue 



-po 



24 



k 2 Wapo ,„ 4 



11520 



fc 4 (20) 



with k 2 = q 2 + p 2 . The term oc fc 4 is always stabilizing, 
but the homogeneous basic state becomes unstable with 
respect to density modulations for a positive coefficient of 
k 2 leading to the corresponding critical filament density 



1 



Pd 



(21) 



providing that a > — holds. The corresponding eigenvec- 
tor is Ed = (1,0, 0) T and the dispersion is shown as the 
dotted lines in Fig. |TJi) and b) for subcritical and super- 
critical values respectively. 

Of the remaining two eigenvalues, 07 and as, only the 
first may become positive in a finite range of k as indicated 
by the solid lines in Fig. ^ It describes an instability with 
respect to orientational fluctuations with the dispersion 
relation 



apo 



768 



apok' 



(22) 



while ag, shown as the dashed lines in Fig. ^ is always 
dampened and related to diffusion of orientational modes. 
The two corresponding eigenvectors are 






E s = I q 



(23) 



It should be mentioned that the eigenvalues aa(k), 
ai(k), as{k) depend only on even powers of the wave num- 
ber modulus reflecting the rotational symmetry of the ba- 
sic state (po, t). The restabilizing fc 4 -term in Eq. I)22ll was 
missing in Ref. |26| . as pointed out recently |27j . However, 
only the interplay between the k 2 and fc 4 contribution as 
in Eq. Ill'2|l allows the identification of a finite wavenum- 
ber instability. The extremal condition together with the 
neutral stability condition 



dai 
~dk 



, aj(k = k c 







(24) 




Fig. 1. The real parts X(k) = Ke[a(k)] of the eigenvalue ad 
corresponding to the instability with respect to density fluctu- 
ations (dotted line) as well as the one with respect to orien- 
tational fluctuations, ai (solid line), are shown as a function 
of the wave number k for D r = 0.15. The third eigenvalue as 
is always dampened as depicted by the dashed line. In part a) 
we have chosen a — 21 and therefore p c < po < Pd leading to 
an orientational instability, whereas in part b) one has a = 26 
and pd < p c < po, hence both modes ad and ai have positive 
real parts in a finite wave number range (see also Fig. El - 
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Fig. 2. The critical densities p c (solid line) and pd (dashed 
line) for the orientational and density instability respectively 
are shown as a function of a and for D r = 0.15. The dash- 
dotted horizontal line is the critical density piN = §7r, above 
which the isotropic-nematic transition due to excluded volume 
interactions takes place. In the range referred to as N one has a 
pure homogeneous transition to nematic order, in range O one 
has a motor driven spatially periodic orientational order and 
in range D one has modulations of the filament density. The 
region S denotes the parameter range where the homogeneous 
solution is stable. In the following sections of the paper we 
investigate the nonlinear behavior in the O-region. 



allow the determination of the critical filament density p c 
at the critical wavenumber k c , above which the orienta- 
tional instability takes place: 



84 
5LL 



K. = 4 



/12A 



1/4 



(25) 
(26) 



Fig. n displays the dispersions, i.e. the wavenumber de- 
pendent growth rates, ad(k), ai(k) and as(k) as dotted, 
solid and dashed lines respectively. In Fig. the criti- 
cal density pd (dashed line) for an instability with respect 
to inhomogeneous density fluctuations and p c (solid line) 
with respect to inhomogeneous orientational fluctuations 
t(r, t) are shown as a function of motor activity a. The 
dash-dotted line describes the critical density p/jv <|12|1 
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above which the homogeneous isotropic-nematic transi- 
tion takes place. On the left side of the dotted line, orien- 
tational fluctuations have lowest threshold, while on the 
right side density fluctuations become unstable at first. 
For increasing values of the rotational diffusion coefficient 
D r , which determines how negative the dispersions 07 (fc) 
and crs(k) start from k = 0, the solid line is shifted up- 
wards in Fig. |3 decreasing the a-range wherein the orien- 
tational instability has lowest threshold with p c < pin- 

Beyond some critical value of (3, the /3-contribution to 
the active current (|8I9|I induces an oscillatory bifurcation 
from the homogeneous filament distribution as described 
in a forthcoming work, where also the nonlinear behavior 
of the oscillatory patterns is discussed . The contribu- 
tion of the rotational current H8I10|I to Eqs. (|f 5JI includes 
only nonlinear terms and therefore does not influence the 
thresholds of the instabilities. Its influence on the non- 
linear behavior of stationary periodic patterns is among 
other things discussed in the next section. 



4 Nonlinear analysis of the orientational 
instability 

The amplitude of the linear solution in Eq. (jf 9J1 is limited 
by the terms in Eqs. H15(l that are nonlinear with respect to 
p and t. Here we investigate the weakly nonlinear behav- 
ior beyond the orientational instability, i.e. in the region 
referred to as O in Fig. |3 where p c < p < pd holds. This 
can be done numerically as exemplified in the next section 
and even analytically if po is not immoderately beyond p c 
defined in (|25[1 : for a supercritical bifurcation, the ampli- 
tude of the mode initially growing with cri{k) is small im- 
mediately above threshold and may be determined in this 
range by a perturbative analysis, the so-called amplitude 
equation [2EJ, whose derivation is sketched in appendix IDl 

The generic form of amplitude equations depends on 
the preferred pattern beyond a stationary supercritical 
bifurcation. In two spatial dimensions the pattern is ei- 
ther spatially periodic in one direction (stripes), in two 
(squares/rectangles) or in three directions (hexagonal pat- 
terns) |28j . Numerical simulations of the basic equations 
<|15|) indicate, as described in Sec. [3 that stripe and square 
patterns are favored immediately above the threshold of 
the orientational instability. Moreover, hexagonal struc- 
tures are not driven in this system due to the overall up- 
down symmetry of Eqs. (| 1 f>|) with respect to t, c.f. ap- 
pendix El as well. 

Square patterns at threshold may be described analy- 
tically by a superposition of two linear modes (i.e. stripes) 
with orthogonal wave numbers which can be chosen with- 
out restriction as ki = (fc c ,0) and k2 = (0, k c ) with k c 
given by Eq. (j26|) leading to the ansatz 



less control parameter 



Wl 




e ik c x + 



o ik oV 



c.c. 



(27) 



where c.c. means the complex conjugate. In a small neigh- 
borhood of the threshold p c , measured by the dimension- 



Po 



(28) 



the time-dependent amplitudes X(t) and Y(t) of the square 
pattern are determined by two generic coupled nonlinear 
equations, the so-called amplitude equations |32II28| 

T 8 t X =sX- 9l \X\ 2 X - g 2 \Y\ 2 X , (29a) 
T Q d t Y = sY- 9l \Y\ 2 Y - g 2 \X\ 2 Y . (29b) 

The coefficients To, g\ and gi are derived in appendix IDl 
The specific system under consideration enters into the 
description only in these coefficients via the parameters 
of the basic equations (|I5fl . namely a, D r and 7, while 
Eqs. i|29|) are universal for a whole symmetry class of pat- 
tern forming systems. 

Apart from the trivial solution Xq = Yq = 0, reflecting 
the stable homogeneous system, the coupled amplitude 
equations (|29|l have also stationary finite amplitude solu- 
tions. These are at first 



X Q = ± 



X = , Yn 



Y = 



(30a) 
(30b) 



which correspond according to Eq. (|27|l to stripes peri- 
odic either in x- or in y-direction. Secondly, there is the 
stationary solution of equal amplitudes 



X = Y Q = ± 



9i + 92 



(31) 



which constitutes a square pattern. In real space, this 
square pattern in terms of the components of the vec- 
tor field t(r, t) resembles the structures which are called 
asters and found in numerous experiments jHJ as will be- 
come clear from the simulation pictures in Sec. El 

As summarized in appendix [El by a linear stability 
analysis of the (nonlinear) stationary solutions given by 
Eqs. I|30|) and (|31|) one finds stable stripes as the preferred 
solution in the range gi > g\ > of the nonlinear coeffi- 
cients. In the parameter range I32I < 5i the square pattern 
is preferred |32) . 

These criteria for g\ and g2 may be translated accord- 
ing to their parameter dependence into the D r -a plane, as 
shown in Fig. [3] with the nonlinear contributions of the ro- 
tational current neglected at first, i.e. 7 = 0. The analytic 
calculations presented here are valid below the solid line 
in Fig. |21 since above the density instability takes place. 
The dotted line corresponds to the condition g\ = g2 > 
which separates the range of stable squares from stable 
stripe solutions. Along the dash-dotted line in Fig. the 
bifurcation from the homogeneous basic state to the stripe 
pattern changes its behavior from a supercritical (below) 
to a subcritical one (above), so the expansion method is 
no more effective. 

Nonlinear effects of the rotational diffusion and the 
active rotational current are described by the last term 
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Fig. 3. The stability regions of stripes and asters are shown 
as calculated by the amplitude expansion method. The two 
critical densities p c and pd coincide along the solid line and 
beyond it instabilities with respect to density modulations are 
preferred not included in our present nonlinear analysis. The 
dotted line is given by gi = g^ and separates the range of 
stable square patterns (asters) from the range of stable stripe 
patterns. Along the dash-dotted line one has g\ — and the 
bifurcation to stripes changes from supercritical (below) to a 
subcritical one (beyond). Between the dashed line, which is 
determined by g 1 — — <?2, and the dash-dotted line asters can 
still exist but are unstable while the amplitudes of stripes can- 
not be determined by our lowest order expansion. Beyond the 
dashed line, also asters bifurcate subcritical. 
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One may complement this outline of the bifurcation 
behavior by a discussion of the two decisive model pa- 
rameters, namely a and D r , and analytical estimates for 
them. Simple models for motor proteins [2TllrjlTT7 26 im- 
ply that the rate a of the translational active transport 
grows linearly with the active motor density p m and with 
the length of the filaments, i.e. a oc p m l- Hence this rate 
can be controlled by the cell in the most effective way by 
the degree of motor activity (i.e. by regulating the ATP 
concentration) as well as on a much larger timescale by 
the density of the motors and the filament length. 

For the rotational and translational diffusion coeffi- 
cients in a dilute solution, calculations taking the hydrody- 
namic interaction into account rT5| propose the analytical 
expressions 



31n(f/6) 
irrjP 



ln(//6) 
2irr]l 



(32) 



In our scaled units this means D' r = -j^D r = 6, lying far 

in the range of squares (asters). 

For semi-dilute solutions one obtains in a similar man- 
ner D r = 6D||/(Z 2 (1 + Crp'f) or 



I 2 6 
D> = —D = - 

r Dn (l + Cp') 2 



(33) 



in scaled units, where c r ~ 1 is a geometry factor from a 
tube model calculation. Since stripes or bundle-like struc- 
tures are stable for D' r < 0.3-0.4, c.f. Figs. Eland 0] one 
needs a rather high (but possible) filament density for such 
a one-dimensional ordering. 

According to these estimates asters are the most likely 
pattern occurring above a stationary bifurcation in dilute 
or semi-dilute two-dimensional motor-filament-systems. 
For bundle-like structures to emerge rather high filament 
densities are needed, which is physically intuitive from the 
overlap nature of all the interactions, namely excluded 
volume and motor-induced filament-filament interaction. 



5 Results of numerical simulations 



Fig. 4. The same diagram as in Fig. [3] but here the nonlin- 
ear contributions of rotational diffusion and active rotational 
currents, c.f. Eq. 1101 . are taken into account with 7 = 0. The 
region of stable stripes considerably broadens while the regions 
of subcritical bifurcations are moved to higher a and lower D r 
values. 



in Eq. i|15bjl . Since they make no contribution to the lin- 
ear operator, they can only change the pattern selection, 
i.e. the stability regions of the nonlinear solutions. Tak- 
ing them into account, the bifurcation behavior from the 
homogeneous basic state is changed as shown in Fig. 0] 
for 7 = a, showing that these effects enlarge the range of 
stable stripe patterns in the D r -a plane. 



Besides the weakly nonlinear analysis described in the pre- 
vious section, the basic equations (|T3|) have been solved 
numerically in order to check the validity range of the per- 
turbation analysis and further explore the solution space. 
For this purpose a Fourier Galerkin pseudo-spectral me- 
thod has been used imposing periodic boundary condi- 
tions on the system. 

Since the validity range of the amplitude expansion 
with respect to the control parameter e is not known a pri- 
ori, in Fig. we compare the amplitude of a stripe solution 
as obtained by a numerical solution of the basic equa- 
tions l(T31) with the analytical results given by Eq. lj3UI) . 
Close to threshold there is nearly perfect agreement be- 
tween both approaches. However the validity range of the 
amplitude equations is actually restricted to a range be- 
low e ~ 0.006 for the parameters used in Fig. [3] Around 
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Fig. 5. A comparison of the amplitude of a stripe solution as 
obtained in simulations (crosses) with the analytical calcula- 
tion Xq — ±(e/gi) ' as given by Eq. igOJ is shown. The a gr ce- 
ment remains well for higher values of e, but in the range be- 
tween e ~ 0.006—0.007 a secondary bifurcation takes place ren- 
dering the high-e solutions unstable. Parameters are D r = 0.5, 
a = 21 and 7 = 0. 



this value a secondary instability takes place, which is not 
taken into account in the perturbation expansion. The nu- 
merical simulations show that beyond this secondary in- 
stability a pronounced accumulation of the filaments to 
densities even higher than pi^ appear accompanied with 
high alternating orientations. These solutions are numer- 
ically stable but nevertheless in an invalid range of our 
model since the nematic order parameter, c.f. Eqs. l(TT)l . 
has been neglected in the moment expansion which is not 
legitimate anymore - so we reject showing pictures here. 

In addition to the validity range with respect to e, one 
may also confirm numerically the stability of the nonlinear 
solutions as predicted in Fig. by the weakly nonlinear 
analysis. As an example, we start with a square solution 
as shown in Fig. a t the point a = 21 and D r = 0.15 
in parameter space (and e = 5 • 10~ 4 ) belonging according 
to Fig. |3 to the region of stable stripe patterns. After a 
slight perturbation, the simulated temporal evolution in 
Fig. Et) shows that only one of the initially equal ampli- 
tudes remains finite in the long time limit leading to the 
predicted stationary stripe pattern displayed in Fig. Et>). 
By several numerical runs we confirmed the analytically 
predicted stability diagrams in Figs. |3] and 

In the parameter range of stable asters, a vectorplot 
of the orientation field superposed on the color coded fil- 
ament density is presented in Fig. The filament orien- 
tation is indicated by arrows, the length being a measure 
of the degree of orientation. The density is high in the 
bright regions and low in the dark ones. At the right hand 
side in Fig. Q one can spot an aster with arrows point- 
ing radially from a center with lowered filament density 
opposing to an inverse aster top left with arrows point- 
ing radially into the center (we remind the reader that 
periodic boundary conditions are imposed). The centers 
of the asters have lowered filament densities which can be 
explained by the nonlinear analysis: in the derivation of 




Fig. 6. A scenario is shown where a repetitive structure of 
asters as displayed in (a) is slightly disturbed at a point in 
parameter space where stripes are preferred, namely at a — 21 
and D r — 0.15 ( and e = 5 ■ 10 -4 ). In part (c) the solid 
lines correspond to the analytical values of the amplitudes 
X s = Y s = 0.0297 for squares and X r = 0.0467 for stripes 
respectively while the dotted lines show how the amplitudes X 
and Y of the polar orientation components t x and t y evolve in 
time from the unstable square pattern to the stationary stripe 
solution in (b). 




Fig. 7. A simulation of Eqs. 115H in the range of stable station- 
ary asters is shown as a superposed plot of the orientation field 
(arrows) and the filament density (dark color coding low den- 
sity, light color high density). Parameters are a = 21, D r = 0.5, 
e = 5 ■ 10" 5 . 



the amplitude equations in appendix [D] one can clearly 
see that the growing amplitudes of the orientation modu- 
lations exite higher density modes - which then limit the 
orientation amplitudes to render the system stable - and 
therefore in the center of an aster, where the orientation 
vanishes, there is no need for a high density. Both density 
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and the degree of orientation reach their maximum in be- 
tween the asters and two saddle-like structures building 
up a square with the two opposing aster centers complete 
the repetitive structure of the pattern at threshold in the 
motor-filament-system. 

6 Summary and conclusions 

For a model of interacting biopolymers and motor pro- 
teins, a nonlinear competition between the formation of 
stripes and square patterns has been described for the 
first time. The model is based on a nonlinear (and nonlo- 
cal) Smoluchowski equation for the distribution function 
of rigid rods extended by active currents to account for 
the motor-mediated relative motion of the filaments. This 
equation can be approximated by a local equation of mo- 
tion under the assumption that the spatial variations are 
small on the length scale of the filaments, which allows 
a gradient expansion |26| . However, this expansion has 
to be continued up to fourth order in the gradients as 
pointed out recently in order to reasonably evaluate the 
linear stability of the homogeneous filament distribution 
|27| . This expansion as well as the linear stability of the 
basic state are described comprehensively. The homoge- 
neous state may become unstable with respect to density, 
orientational and nematic instabilities, but in different pa- 
rameter ranges as shown in Fig. |5] The two crucial model 
parameters are the rate of active translational transport 
and the rotational diffusion coefficient which both can be 
regulated by the cell by creating filaments, motor proteins 
and the fuel ATP as discussed and estimated at the end 
of Sec. H 

Our subsequent nonlinear analysis focuses on the ori- 
entational instability, which may be of high biological rel- 
evance for aster-like structures as in the mitotic spindle 
which have also be seen in in vitro systems recently 
This analysis is to a large extent analytically in terms of an 
amplitude expansion technique yielding two coupled equa- 
tions, namely Eqs. (|29() . for the amplitudes of the spatial 
modulation of the orientational field in x- and y-direction. 
These two coupled equations determine the amplitudes 
of square (asters) and stripe patterns and they provide 
also the stability boundary between both solutions near 
threshold. Simulations with a spectral code and periodic 
boundary conditions confirmed these analytical results. 

Estimates for the two model parameters a and D r 
of the basic equations (|15p. c.f. end of Sec. 0] suggest 
that asters are the most likely pattern occurring above a 
stationary bifurcation in dilute or semi-dilute two-dimen- 
sional motor-filament -systems. Stripes may be expected 
only for rather high filament densities which can be traced 
back to the overlap nature of both present interactions, 
namely excluded volume and motor-induced filament-fila- 
ment interaction. Once again both the type of the pattern, 
i.e. asters or stripes, as well as its wavelength could be gov- 
erned by the cell by creating filaments, motor proteins and 
ATP. 

Experiments as described in Refs. 7,8,5] show that in 
a model system comprised of microtubules and a single 



type of motor in a confined geometry asters may arise but 
are unstable to a vortex with motors rotating around the 
center, whereas in unconfincd geometries, with increas- 
ing motor concentration, the arising patterns are vortices, 
mixtures of asters and vortices, asters and finally bun- 
dles of microtubules for high motor densities. This is in 
agreement with our calculations visualized in Fig. as 
for moderate values of D r the asters in our model become 
unstable against stripe formation for increasing a, which 
is proportional to the (homogeneous) density of motors. 
Since we only considered periodic boundary conditions, 
we did not find vortices, which nevertheless may be possi- 
ble if a confined geometry is simulated. In contrast to this 
good agreement with microtubule-motor systems, recent 
experiments on actin-myosin systems showed that the 
models at hand can not be directly used in that case. But 
modifying the present model to account for ATP deple- 
tion and stable cross-linking of the filaments by inactive 
myosin complexes, it is even able to explain the different 
aster formation mechanism in actin |11I35| . 

The first attempt of describing patterns in microtu- 
bules and motors in Ref. [Hj was done by a molecular 
dynamics simulation where asters and vortices could be 
reproduced. In a subsequent publication |12| . the case of 
two types of motors walking on the filaments in opposite 
directions was investigated and a pattern referred to as 
a network of interconnected poles was obtained which is 
very similar to the pattern near threshold displayed in 
Fig. [7] A first macroscopic description of the problem was 
given in Ref. jI4) , where coarse-grained equations derived 
from a spin- like model displayed stripe patterns. Finally, 
in Ref. ^B] equations for the motor density and a phe- 
nomenological vector field describing the polar orientation 
were proposed, which reproduced asters and, in a confined 
geometry, vortices, but were not able to describe stripe 
patterns. This can be traced back to the omission of an 
explicit equation for the filament density. More severely, 
as will be worked out in Ref. |35|. in the model presented 
here linear couplings between the density and the orienta- 
tion field via the active transport terms proportional to (3 
lead to an oscillatory instability and propagating modes 
through the motor-filament-system, which cannot be de- 
scribed by Ref. JS] but should be highly relevant for vital 
processes like cell spreading on a substrate as measured 
for instance in Ref. |36| . 

Some experiments indicate that the spatially inhomo- 
geneous distribution of the motor density may become 
a further relevant degree of freedom and therefore 
should be taken into account for various purposes in fu- 
ture models. This degree of freedom in combination with 
stable crosslinking may also give rise to a different insta- 
bility in actomyosin 11,35]. From the theoretical point of 
view it would be interesting to extend the moment expan- 
sion scheme one step further yielding an equation for the 
nematic order parameter tensor. This would allow to take 
also the effects of nematic order into account and to study 
its interplay with the polar orientation in the presence of 
motors, which may lead to an interesting model system of 
a new symmetry class. 
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A Symmetries of the motor-mediated 
velocities 

The explicit form of the motor-mediated translational and 
angular velocities, namely Eqs. © and I jlOjl . can be ob- 
tained by writing down the simplest terms fulfilling the 
conservation laws and symmetries of the system. 

If one considers an interacting filament pair in the 
absence of external forces and torques, both momentum 
and angular momentum of the pair have to be conserved. 
Hence the motor-mediated translational and angular ve- 
locities, v = r — r' and u> = u — u' , have to be odd under 
the transformation (r,u;r',u') — ► (r',u';r, u). Transla- 
tional invariance leads to /(r, r') = /(r — r') for / = v, u). 
Together with the rotational invariance, v has to be odd 
and u> even under the transformation (r — r';u, u') — > 
(V — r; — u, — u'). The simplest terms fulfilling the above 
conditions for v are proportional to r — r' or u u' and 
for u) proportional touxu'. One can thus write down 

/ , m ar-r'l + u-u' j3 u' - u 
V (-^ U ')= 2^T^ + 2^7l < (34) 



a>(u,u') = 7(11 • u') 



U X u 

|u x u' 



(35) 



The term proportional to a is made orientation-dependent 
(it could be generalized by allowing two different prefac- 
tors) and the plus sign favours parallel alignment. In the 
term proportional to 7, spatial dependence has been ne- 
glected and the prefactor uu' models the tendency of mo- 
tors to bind on two filaments which share an angle smaller 
than -|. The common factor |u x u'| _1 is just a normal- 
ization to the interaction volume. 



by the four dimensional integral 

L/2 L/2 

V ex (r,u) = J du 'J dr ' J d C J d V \uxu'\ 

-L/2 -L/2 

5 (r - r' + uC + u'77) #(r', u') , (37) 

which may be approximated by a Taylor series expansion 
of the (5-function with respect to 77 and £. 2 Performing the 
77 and £ integration one finally obtains 



K x (r,u)=L 2 y"du'|uxu'|- 

l + g{(u.9 r ) 2 + (u'.a r ) 2 } 



!?(r,u'),(38) 



where terms up to second order in the spatial derivatives 
have been taken into account. The prefactor L 2 reflects 
the two-dimensional excluded volume. 

By a similar procedure, the contribution to the active 
current proportional to a can be evaluated as 

xL 3 



J? = -^t di J ^(r,u)(l + u-u') 



24 

Ui(u ■ d r ) 

Win'- dr) 



<f(r,u'). 

(39) 



One should note that due to the f, ^-integrations, odd 
powers of d r vanish in both expressions, which is the math- 
ematical reason underlying the rotational symmetry of 
Eqs. US . 



B The gradient expansion of the interaction 
integrals 

The excluded volume interactions, c.f. Eq. ©, as well as 
the motor induced rod-rod interaction, c.f. Eqs. (jSJ, are 
defined by overlap integrals. Hence the equation of mo- 
tion {2J f° r the probability distribution function ^(r, u) is 
nonlocal and its solution is exceedingly difficult. Assum- 
ing that spatial variations are small on the length scale of 
a filament, in order to deal with a local equation one can 
perform a systematic expansion of the integrals with re- 
spect to gradients of the probability distribution function, 
as described in this appendix. 

The interaction kernel given by Eq. (JJJ) can be ex- 
pressed in terms of Straley coordinates |33j , which in two 
dimensions are defined by 

r — r' = \iQ + u'77 1 (36) 

with the parameter constraint < £, 77 < L/2 (L the 

filament length) and the Jacobian |u x u'|. Thus the ex- 
cluded volume interaction, namely Eq. ©, is determined 



C The moment expansion method 

The method for deriving the coupled Eqs. (jl5f) for the den- 
sity and the orientation field of the filaments from the un- 
derlying Smoluchowski-equation, Eq. , is similar to the 
calculations presented in Ref. |SH1 , where it has been used 
to obtain equations describing the initial stage of spinodal 
decomposition during the isotropic-nematic transition of a 
three-dimensional hard-rod fluid. In both cases, the start- 
ing point is an approximate description of the probability 
distribution function "^(r, u, t) by its moments, but in con- 
trast to a usual liquid crystal as considered in Ref. |38| . 
now the rods are polar with respect to the motor action. 
Thus the ±u-symmetry is broken and the first moment, 
namely the filament orientation field t, does not vanish 
here and cannot be omitted. Hence one has to approxi- 
mate !^(r, u, t) by its first two moments, 

iP(r, u, t) g ^ {p(r, t ) + 2u • t(r, t )) , (40) 

2 Equivalently, after performing the r'-integration the 
shifted distribution function & (r + u£ + u'77) may be ex- 
panded. 
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with the filament density p(r, t) and the orientation field 
t(r,i), c.f. Eqs. fTTJi . 

The validity of the representation in Eq. (|40|) can be 
seen immediately by using it to evaluate the first mo- 
ments which correctly yields J du ^(r, u) = p(r) and 
J (iuu!f(r, u) = t(r). Inserting the currents as defined by 
Eqs. J2J and JSJ as well as applying the gradient expan- 
sion given in appendix^! then an integration of Eq. (J2J by 
J du and J duu yields two evolution equations for p(r,t) 
and t(r, t) respectively, as given by Eqs. fT5|l . 

The remaining integrals involved are merely orienta- 
tional averages. If one defines 



d6 
2^ 



A{6) , (41) 



where 9 parameterizes the unit vector u in two dimensions, 
the following formulas are useful and easily proven: All 
mean values depending on odd powers of u vanish as well 
as the mean values of any product of |u x u'| and odd 
powers of u or u', since |u x u'| = y/l — (u- u') 2 contains 
only even powers of u, u'. For the even powers one gets 



(u a up) u = -5 a/ 3 , 



(42) 



{UaU^U^U^u = - ((5 Q/3 <5 M „ + 8 al j,bf}v + fiavSflfj.) (43) 



and 



{U a Uf 3 U ll U p 



v U a U T ) u = — ^2Sij5klSmn , (44) 



where V means all the permutations of {a, /?, /i, v, a, t} 
generating different index combinations of the Kronecker 
delta product. Furthermore one needs 

<|uxu'|) u ,= / ^\ 8 in(9-6')\ = - (45) 

J Z7T 7T 



and 



2 

(|u x u.'\u' a u'p) u , = -— [u a u - 28 a0 ) . (46) 



The operator of rotational diffusion in two dimensions as 
given by Eq. (J5J can be expressed as 

[K]i = [u x d u ]i = 8 i3 (uiu' 2 - u 2 u'x) (47) 

and for the rotational terms, the following partial integra- 
tion formula 

(A(u)K [£(u)]) u = —(H [A(u)] B(u)) u (48) 

is crucial to simplify calculations. 

As an example, we calculate a rather simple term, 
namely the excluded volume contribution to the equation 
of motion for the filament density p, c.f. H15a(l . which is 



proportional to Dm and second order in spatial derivatives. 
So we have to deal with 

L>H Jdu u iUj di (&(u)djL 2 Jdu'\u x u'|^(u')^ , (49) 

where spatial coordinates have been suppressed. Substi- 
tuting the representation Eq. (|40|l twice and using the 
nomenclature introduced in Eq. Ij41(l , we have to calculate 

(u^di ([p + 2u k t k ] dj(\u x u'| \p+ 2 U ;t ; ])) u ) u ,. (50) 

The contributions u^tk, u\ti vanish because of their odd- 
ness. Thus we are left with 

(uiUjdi (pdj{\u X u'|) u /p)) u (51) 

which can be simplified to 

-(uiU^ndi {pd ]P ) = -V • (pVp) (52) 

7T 7T 

by using Eqs. |@2J| and l|43|) . 

D Derivation of the amplitude equations 

Here we describe the scheme for the derivation of the two 
coupled amplitude equations (|29Jl from the three under- 
lying nonlinear equations (|15fl . First of all one assumes 
small values for the amplitudes X and Y of the spatially 
periodic deviations from the homogeneous basic state po 
and t = 0. At threshold, i.e. at po = p c and for k c as cal- 
culated in Eqs. (|25f) and Q26p. these deviations are either 
periodic in x- or in y-direction as described by the ansatz 
in Eq. |27fr. 

Similar to Sec.|3| the nonlinear equations l|15|l may be 

rewritten in terms of the deviations w = (p,t Xl t y ) from 
the basic state w = (po, 0, 0) as follows: 



with 



<9 t w 



N(p,t) 



£ w + N(p,t) 

^ p ( P ,t) N 
A4(p,t) 

MM j 



(53) 
(54) 



The linear operator is defined by Eq. I|18|) and the non- 
linear operator N includes all the nonlinear terms with 
respect to p and t on the right hand sides of Eqs. I|15|) . 

Naturally, as the small expansion parameter the rela- 
tive distance from the threshold, 



Po 



(55) 



is chosen. Close to threshold the dynamics of the linear 
solution in Eq. I|19f) is slow and accordingly a slow time 
scale 



T = et 



(56) 



is introduced allowing the time derivatives in Eqs. I|15|) to 
be replaced by 



(57) 
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The solution w is expanded with respect to powers of e 1 ^ 2 
w = e 1/2 wi + ew 2 + e 3/2 w 3 + . . . , (58) 
with wi as in Eq. (|27Jl . as are the nonlinearities 

N = eN 1 +£ 3 / 2 N 2 + ... . (59) 



Sorting the contributions to Eq. i|53|) with respect to pow- 
ers of £, one ends up with the following hierarchy of equa- 
tions 



,1/2 

e 

.3/2 



£oWi = , (60a) 
C w 2 = -Af p (ti)e p , (60b) 

£qw 3 = 3 T wi - £ 2 wi -} A/i(p2jti)ei , (60c) 



that has to be solved successively. We have introduced 



(1,0,0), e x = (0,1,0), e y = (0,0,1) and 



Hi = 



o 

V o c 





r (2) r (2) 

/ -22 / -23 

(2) £ (2) 



(61) 



■ VI 



-33 



with 



(2) 



Ml 

Po 

—C {2) 
Po 
1 



— — — ) _\ 

2tt 24, 



-22 



(2) 
23 



a 
96 
a 
"48 



(A + 2dl 

d x d, 



19 a 
11520' 

a 

46080 



Zi 2 



(llZ\ 2 + 64A9 2 



720 a 



(62) 



The remaining two matrix elements £32'' and follow 
from £33 and £22 by permuting d x and d y . 

The equation in C(e 1//2 ) is just the linear eigenvalue 
problem already discussed as Eq. (fTp|l in Sec. , i.e. it is 
solved by pi = and 



'(2) 



t lx = X(T)e lk " x + c.c. , ti v = Y(T)e lk <= v + c.c. (63) 

describing the orientational fluctuations. 

At the next order 0(e) , the nonlinearity is only present 
in the density equation, while £0 acting on the t-subspace 
is nonsingular, leading to t 2 = 0. The overall up-down 
symmetry of Eqs. (|15f) with respect to t prohibits hexag- 
onal structures as can be seen in this order in s, where 
for hexagons to be driven a quadratic contribution in the 
equations that are linearly unstable, i.e. in the equations 
for the orientation fields t x and t y would be needed |31| . 
Inserting ti in J\f p (t{) yields an equation for p 2 , whose 
solution is of the following form 

p 2 (X, Y) = ri X 2 e 2lkx + r 2 Y 2 e 2lqy 

+ r 3 XYe l(kx+qv) + ri XY*e l(kx ' qy) + c.c. (64) 

with n = Ti{a, -D r ,7) for i — 1,..,4. Here one can see 
that the coupling of the orientational field to the density 
is crucial for the physical stability of the system, since the 



density p 2 is responsible for the saturation of the ampli- 
tudes in the equation of the next order e 3 / 2 , c.f. Eq. (|60c|) . 
while t 2 — and thus t can not limit the amplitudes. 

Instead of solving Eq. I|60c|l at the order 0(e 3 / 2 ), one 
can use Fredholm's alternative, which states that for Eq. 
(|60c|) having solutions, there must not exist terms on its 
right hand side that lie in the kernel of Co, i.e. no con- 
tributions proportional to the critical modes exp(ik c x), 
exp(ik c y). Collecting the prefactors of these respective 
modes, one gets the two equations (1291) . with analytical 
but lengthy expressions for to, g\,gi as functions of a, D r 
and 7. 



E Existence and stability of the nonlinear 
stripe and square state 

The range of existence as well as the range of linear sta- 
bility for roll solutions and for the squares can be easily 
investigated in terms of the amplitude equations. Station- 
ary single amplitude solutions as given in Eqs. l|3T))l exist 
beyond threshold, i.e. for e > 0, only when g\ > holds. 
The amplitudes for a stationary square, Eq. I|31(l . follow 
from Eqs. (|29J) assuming equal amplitudes 



\X\ = \Y\ = 



.91 + 92 



(65) 



Obviously, squares exist beyond threshold only in the pa- 
rameter range gi + g% > 0. 



Linear stability. For small perturbations SX and 5Y of 
the stripe and square solutions respectively, one obtains 
by the ansatz X = Xo + SX and Y = Yo + 5Y and lin- 
earizing Eqs. H29|) with respect to 8X and SY two coupled 
equations in both perturbations. Those may be solved by 
the mode ansatz (SX, SY) ~ (SX, 8Y)e at leading to a sec- 
ond order polynomial in a providing two eigenvalues. One 
of them is always negative while the second is either 



9i -92 
9i 



for rolls or 



(is = 2e 



92 - ffi 

gi + g2 



(66) 



(67) 



for squares. 

Thus stripes or squares are stable if a r or a s is neg- 
ative, respectively. Accordingly stripes are the preferred 
solution in the range of the nonlinear coefficients g 2 > 
gi > 0, while in the parameter range \g%\ < g\ the square 
patterns are preferred |32j . Since the nonlinear coefficients 
<?i and gi are functions of the rate of active translational 
transport a and of the rotational diffusion coefficient D r , 
we are able to plot the stability ranges of the patterns as 
depicted in Figs. El HI 
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